A Systematic Review of INGARCH Models for Integer-Valued Time Series

Count time series are widely available in fields such as epidemiology, finance, meteorology, and sports, and thus there is a growing demand for both methodological and application-oriented research on such data. This paper reviews recent developments in integer-valued generalized autoregressive conditional heteroscedasticity (INGARCH) models over the past five years, focusing on data types including unbounded non-negative counts, bounded non-negative counts, Z-valued time series and multivariate counts. For each type of data, our review follows the three main lines of model innovation, methodological development, and expansion of application areas. We attempt to summarize the recent methodological developments of INGARCH models for each data type for the integration of the whole INGARCH modeling field and suggest some potential research topics.


Introduction
This paper reviews the development of modeling and inference for four types of integer time series, including the unbounded N -valued counts, Z-valued time series, bounded N -valued counts and multivariate counts. Firstly, the unbounded N -valued counts, also known as count time series, refer to discrete time series taking values in the range N = {0, 1, 2, · · · }. For example, during an influenza outbreak, the number of new confirmed cases is reported daily, even down to each community. The analysis of such data is one of the fundamental tasks in epidemic forecasting and policy implementation (see Agosto and Giudici [1], Agosto et al. [2] and Giudici et al. [3] among others). Secondly, Z-valued count time series taking values in the range Z = {· · · , −1, 0, 1 · · · } are the appropriate tool to employ when attention is turned to, for example, the changes in athletic performance by the difference between the number of goals scored in each game and that in the previous one. Thirdly, the study of bounded count time series with a range of {0, · · · , n} for a given n ∈ N is also a concern, such as the rigorously recorded data set of water quality in the estuary. Finally, with the development of data acquisition technology and the expansion of storage space, multivariate count time series have emerged in various fields, and the related research is also developing.
For the above four types of integer time series, this review focuses on the the relevant research fields of the integer-valued generalized autoregressive conditional heteroscedasticity (INGARCH) models, which assume that the observations follow a discrete distribution and are conditioned on an accompanying intensity process that drives the dynamics, which are classified as observation-driven models. For instance, given the intensity process, the observations of a count time series follow a Poisson distribution and the intensity process is a linear combination of its lagged values and lagged observations (Fokianos et al. [4]). The diversity of count time series urges the INGARCH-type models to evolve to a broader domain. This is accompanied by the development of many subdivision directions, including the selection and even the creation of suitable conditional distributions, the exploration of nonlinear dynamic structures, the proof of stationarity and ergodicity, the relaxation of restrictions, the updating of methodologies, and so on. The purpose of this paper is to review recent developments, including the above subdivision directions, in INGARCH models over the past five years, as earlier work has been summarized. For the generally methodological development of count time series, not only INGARCH-type models, please refer to Weiß [5] and Davis et al. [6].
The INGARCH-type models for count time series are the most well-developed of these four types of integer-valued time series modeling, and they are reviewed in Section 2. On the one hand, to model some specific features of count time series, researchers have innovated the conditional distribution assumption. Indeed, Gorgi [7] and Qian and Zhu [8] focus on heavy-tailed count time series, while Silva and Souza [9] is more concerned with the generality of distribution and Souza et al. [10] further extends Silva and Souza [9] to the time-varying version. Accordingly, they all propose new conditional distribution assumptions that can model the features of interest. On the other hand, there are some new developments in the dynamic structure of INGARCH. For example, Weiß et al. [11] introduces the softplus function in the linear dynamic structure to implement the possibility of negative ACF; similar to Souza et al. [10], Roy [12] also focuses on the time-varying feature, but proceeds from the construction of a semi-parametric dynamic structure based on a Bayesian framework; some advances are made in the study of nonlinear dynamic structures with the help of threshold structures and the hysteresis model (Liu et al. [13], Chen and Khamthong [14], Liu et al. [15] and Chen et al. [16]). Moreover, when it comes to proving the existence and uniqueness of stationary distributions, which has been a challenging work in this field, the common approaches recently are based on approximation techniques, weak dependence or theoretical frameworks using the Feller property, e-chain and Lyapunov methods. In contrast, the absolute regularity for nonlinear GARCH and INGARCH models under a mild assumption are considered by Doukhan and Neumann [17], and a series of methodological studies are specifically reviewed in Section 2.
In contrast, the INGARCH-type models for the other three integer-valued time series are all still areas of ongoing research. The development patterns of all three areas are similar to that of the above-mentioned INGARCH-type models for count time series, including innovation in conditional distribution assumptions, flexibility in dynamic structure and methodological establishment, etc., but the challenges faced by each are different. For Z-valued time series, the difficulty of finding a suitable conditional distribution can be overcome, since it can be handled directly with the help of a sequence of binary or ternary random variables (Hu [18], Hu and Andrews [19] and Xu and Zhu [20]). However, the INGARCH-type models for Z-valued time series require more complex and subtle dynamic structures to implement, and the corresponding theoretical proofs are of higher difficulty. For bounded counts with range {0, · · · , n} for given n ∈ N , two research frameworks are currently divided into those based on vector form (Fokianos and Truquet [21]) and those based on scalar form (Weiß and Jahn [22], Chen et al. [23], Chen et al. [24] among others). The challenge encountered in the former is the computational stress of the estimation of matrix-type unknown parameters, while the difficulty in the latter lies in the rarity of suitable distributions, which are required to be discrete distributions with a range of values, and the accompanying proof of stationarity. Finally, the development of INGARCH-type models for multivariate numerical time series is still in its nascent stage. Lee et al. [25], Kim et al. [26] and Cui and Zhu [27] focus on the INGARCH-type models for bivariate count times series, but, subject to the development of multivariate Poisson distribution, more multivariate models related to these studies have not been extended. Another emerging area is the INGARCH-type models for time-varying network data (Armillotta and Fokianos [28], Armillotta et al. [29], Armillotta and Fokianos [30] and Tao et al. [31]), and many of the attempts that come with the challenge are worthwhile, such as optimization of conditional distributions, nonlinear dynamic structures and time-varying network data with upper bounds on the number of edges.
The rest of the paper attempts to review the INGARCH-type models, which are presented below. Sections 2-5 present the INGARCH-type models for unbounded Nvalued counts, Z-valued time series, bounded N -valued counts and multivariate counts, respectively ( Figure 1). Since there is much content related to unbounded N -valued counts, Section 2 is divided into three subsections to elaborate on the latest advances in models, methodologies and applications, but the other sections are also arranged in this way although they are not split again. Finally, we present some potential research topics from a personal perspective in Section 6.

Count Time Series
We consider count time series in this section, the possible outcomes of which are contained in N . Let {X t } t∈Z be a count process, x 1 , · · · , x T be a finite set of observations, and t be defined as σ-field generated by observations up to and including time t.
The INGARCH(p, q)-type models with p ≥ 1 and q ≥ 0 require two types of assumptions: first, a distributional assumption for X t conditioned on X t−1 , X t−2 , · · · is needed to guarantee the discrete range of X t ; secondly, the conditional mean M t = E(X t |X t−1 , · · · ) is required to be a linear or nonlinear expression in the last p observations and the last q conditional means, thus constructing a dynamic structure. For example, the Poisson INGARCH model assumes that X t , conditioned on X t−1 , X t−2 , · · · , is Poisson distributed with intensity parameters M t and This section begins with a review of recent modeling work on improving these two types of assumptions separately.

Recent Advances in Assumptions of Conditional Distribution
This subsection begins with the progress of research on count time series with outliers. For the INGARCH models, the common approach is to assume a heavy-tailed conditional distribution. Gorgi [7] introduced a heavy-tailed mixture of negative binomial distributions, known as the beta-negative binomial (BNB) distribution, as an distribution of X t conditioned on X t−1 , X t−2 , · · · . See Table 1 for the specific definition of BNB. Relatedly, Qian and Zhu [8] was also concerned with heavy-tailed count time series and thus uses the generalized Conway-Maxwell-Poisson (GCOMP) distribution (in Table 1), which has one more parameter than the Conway-Maxwell-Poisson distribution, but provides a unified framework to handle over-or under-dispersed, zero-inflated, and heavy-tailed count data.

GCOMP (Qian and Zhu [8])
The PMF of GCOMP is is and

MP (Silva and Souza [9])
Let Z be a continuous positive random variable belonging to the exponential family with density function given by where b(·) is a continuous three-times differentiable function and ζ 0 is such that b (ζ 0 ) = 1 and c(·; ·) is a function that maps Then, X belongs to the class of mixed Poisson distributions. If Z follows a gamma distribution with mean 1 and dispersion parameter φ, we find that X follows a negative binomial distribution with parameters µ and φ. Its probability function is given by

Definition 1.
The following model is denoted by GCOMP-INGARCH(p, q), According to the properties of the GCOMP distribution, the approximate conditional expectation and variance are given by and these are fundamentally the keys to the flexibility of the GCOMP-INGARCH(p, q) model. Qian and Zhu [8] also established some properties by assuming that model (2) is approximately stationary. Additionally, Silva and Souza [9] proposes a general class of INGARCH models by introducing the mixed Poisson (MP) distribution proposed by Barreto-Souza and Simas [32]. This not only enriches the distribution types of INGARCH models, but also extends the negative binomial INGARCH process proposed by Zhu [33], and can even further evolves new models such as Poisson inverse Gaussian and Poisson generalized hyperbolic slope processes. Related works include Manaa and Bentarzi [34], Almohaimeed [35], Almohaimeed [36], and Cui and Wang [37], where Manaa and Bentarzi [34] focused on the expansion of time-varying parameters, which extends the time-invariant negative binomial INGARCH(1,1) studied by Zhu [33] to the periodic negative binomial INGARCH(1,1) model. Besides the flexibility of distributions, Silva and Souza [9] construct an expectationmaximization algorithm for estimating the parameters, in particular, the dispersion parameter. This provides a new framework for parameter estimation of the INGARCH-type models, which we believe also leaves a wide scope for further research.
Note that the dispersion parameter φ in Definition 2 is independent of time t, and thus Souza et al. [10] is further extended for this setting. Indeed, the dispersion parameter is made time-varying, φ t , by creating a dynamic structure for it similar to that for the intensity parameter. The resulting model in Definition 3 is called a linear time-varying dispersion INGARCH (tv-DINGARCH) model (nonlinear tv-DINGARCH(p 1 ,p 2 , q 1 , q 2 ) is omitted here). An interesting feature of the linear tv-DINGARCH processes is that to some extent it is more analogous to the original GARCH model than other INGARCH models. In particular, the mean of this model can be constant, while the variance depends on time as in an ordinary GARCH model. This feature is not possible to be accommodated by the standard INGARCH. Hence, they refer to the model as a pure INGARCH process to highlight the degree of association with GARCH that distinguishes it from other INGARCH models.
In addition, zero-or zero-one-inflated INGARCH models are valuable in applications such as insurance. Lee and Kim [38] proposes more general multiple values-inflated INGARCH models in Definition 4. The conditional distribution q(·|η) containing only one unknown parameter is indeed a streamlined setup, but it cannot be ignored that the multiple inflated values will be accompanied by a simultaneous increase in the number of parameters ρ m and the resulting complexity can hardly be offset. From our personal point of view, a possible reason for such a dilemma is that the idea of a multi-inflated model constructed by superimposing indicator functions remains fundamentally indistinguishable from that of the zero-inflated model.
An alternative approach to model generalization is to weaken the assumptions of the structural model allowing for more generalized link functions, and Wechsung and Neumann [39] make some contributions to the estimation of linkage functions. Wechsung and Neumann [39] considered a nonparametric version of the INGARCH(1,1) model, where the link function in the recursion for the variances was not specified by finite-dimensional parameters. This work completed the asymptotic analysis based on the mixed property, benefiting from the application of powerful exponential tail bounds in connection with a chaining argument. This work is highly theoretical, and shows in principle that a sufficiently regular link function of an INGARCH(1,1) process with hidden intensities can be estimated using a nonparametric least squares estimator, where the estimate is chosen from a truly nonparametric class of candidate functions. Further, the consistency rate of the estimator was shown to be nearly optimal. For practical purposes, Wechsung and Neumann [39] also reported an approximate version of the theoretical estimator.
The assumption about the distribution has also been approached from another perspective. Instead of specifying the conditional distribution, Aknouche and Demmouche [40] presents a double mixed Poisson autoregression whose conditional distribution is a superposition of two mixtures of Poisson distributions. It is more flexible compared to Poisson mixtures at the cost of just a few additional parameters. Further, Doukhan et al. [41] considered the mixture of nonlinear Poisson autoregressions, and Mao et al. [42] proposed a more general mixture INGARCH model, which includes s negative binomial mixture INGARCH of Diop et al. [43] and generalized Poisson mixture INGARCH models.

Recent Advances in Dynamic Structures
This subsection reviews the latest advances in dynamic structures of INGARCH. Many recent works have made attempts in this area, but due to space limitations, we mainly focus on the essential dynamic structural innovations, some of which will not be mentioned one by one; for example, Souza et al. [10] also considered the inclusion of covariate/exogenous time series in tv-DINGARCH model, which enhances the applicability.
First, the INGARCH model with linear dynamic structure has been limited in its setting. Specifically, since the mean of a count random variable is a positive real number, the constraints a 0 > 0 and a 1 , · · · , a p , b 1 , · · · , b q ≥ 0 have to both hold in (1). This severely dampens the possibility of a negative ACF. The existing log-linear INGARCH model is an implementable solution to achieve negative ACF values, but at the cost of losing the linear conditional mean and the ARMA-like autocorrelation structure. Resolving the dilemma between the linear dynamic structure and a wider range of achievable ACF values is one of the key steps in the further development of the INGARCH-type models. The contribution of Weiß et al. [11] is to resolve this dilemma with the help of the softplus function. The definition of softplus INGARCH model is given as follows.
The softplus function sp c (x) avoids the drawback of not being differentiable in zero, while being approximately linear for x > 0. The excellent properties of the softplus functions play a key role as they coincide with the breakthrough of the dilemma of INGARCH mentioned above.
What follows is a progression of time-varying features in INGARCH, which has been mentioned in the previous introduction of linear tv-DINGARCH processes in Definition 3, but Souza et al. [10] approached it from the distribution. Roy [12] proceeded from the construction of a semi-parametric dynamic structure, and proposed a time-varying autoregressive models for count time-series based on a Bayesian framework (Definition 7). This is the first attempt to model possibly autoregressive count time series with time-varying coefficients, and the success of this attempt can be attributed in part to the Bayesian framework.

Definition 7 (time-varying Bayesian INGARCH model).
If the conditional distribution for X t given t−1 is Poi(λ t ) and where the hierarchical prior on unknown functions µ(·), a i (·) and b i (·) are based on B-spline bases, {X t } is said to follow from the time-varying Bayesian INGARCH model.
Take (1) with p = q = 1 for example, the essential role of the original dynamic structure is to drive the time-varying nature of the conditional mean M t , but where the parameters ω, α 1 and β 1 do not vary with time. This means that the dynamic structure is allowed to take into account the heterogeneity changes at t − 1 moments and the effects over the entire time period reflected by the time-constant parameters. In contrast, the semi-parametric time-varying dynamic structure (9) corresponds to further strengthening the time-varying property while weakening the average effect over the whole period. This modeling idea is better suited for rapidly changing count time series. In conclusion, Roy [12] contributes a new idea to the study of INGARCH from a Bayesian viewpoint. Its framework for the study of time-varying Bayesian INGARCH models is adapted to general non-stationary time series.
Next are some advances in the study of nonlinear dynamic structures with the help of threshold structures. The classical two-regime threshold autoregressive model allows for many properties associated with nonlinearity, and thus, both Liu et al. [13] and Chen and Khamthong [14] draw on this classical model to implement a nonlinear study of the dynamic structure of INGARCH, with the difference that the latter is based on the Markov switching approach and introduces covariates. It is worth mentioning that Lee and Hwang [44] proposed a generalized regime-switching INGARCH(1,1) model in a fashion similar to, yet different from, Chen and Khamthong [14], which also employs Markov chains with time-varying dependent transition probabilities. The difference is that Lee and Hwang [44] derives recursive formulas for the conditional probabilities of regimes in Markov chains given past information, starting from the Poisson parameters of the INGARCH(1,1) process.

Definition 8. The Markov switching INGARCHX model with state variables is defined by
where s t follows a first-order Markov chain with the following transition matrix: ) is a non-negative parameter vector in state i. Naturally, changing the form of λ t yields the threshold INGARCHX as follows: where Y t−d is the threshold variable determining the dynamic switching mechanism of the model, d is a delay lag and c is the threshold value.
The segmented dynamic structure may lead to sudden changes in the probability of the INGARCH process, which is a crux that needs to be improved. Li et al. [45] considered a hysteretic process with the hysteresis variable {S t } and the hysteresis zone (r L , r U ], which makes the regime-switching mechanism more flexible. On this basis, Liu et al. [15] improved the self-excited threshold negative binomial autoregression (TNBAR) of Liu et al. [13] and proposed the self-excited hysteretic negative binomial autoregression (SEHNBAR) with the hysteresis variable S t = X t . Definition 9. Let {N t , t ∈ Z} be a sequence of the i.i.d. negative binomial process given in Liu et al. [13], then {X t } is said to follow the SEHNBAR model, if R t is the regime indicator with the hysteresis variable Y t−1 , and (r L , r U ) are boundary parameters of the hysteresis zone satisfying 0 ≤ r L ≤ r U < ∞.
λ t is at the lower regime when X t−1 < r L , while at the upper regime when X t−1 > r U . If X t−1 falls within the hysteresis zone, the regime indicator remains unchanged, which means that the regime indicator at time t is the same with that at t − 1. Formally, the hysteresis model with piecewise linear structure enjoys a more flexible regime-switching mechanism. Another less intuitive but proven advantage is that the lagged observations on which λ t relies are infinitely far away, which is the key difference between hysteresis models and traditional threshold models. Chen et al. [16] also proposed a similar INGARCH model based on the Bayesian framework.
Aknouche and Scotto [46] proposes a multiplicative INGARCH model, which is defined as the product of a unit-mean integer-valued i.i.d. sequence, and an integer-valued dependent process defined as a binomial thinning operation of its own past and of the past of the observed process. This model combines some features of the INGARCH, the autoregressive conditional duration, and the integer autoregression processes, so it can be used to model high overdispersion, persistence, and heavy-tailedness. Furthermore, Weiß and Zhu [47] propose an integer-valued analog of multiplicative error models based on a multiplicative operator, and the resulting models are closely related to the class of INGARCH models.
Last but not least, some migratory research works that have received attention in other fields but are poorly known in the INGARCH model are also of interest. Sim et al. [48] established the overall framework for the study of general-order INGARCH(p, q) models without the restriction p = q = 1. Similarly, the purpose of Tsamtsakiri and Karlis [49] was to select the most appropriate order of INGARCH(p, q) using a trans-dimensional Bayesian approach, and Tian et al. [50] focused on order shrinkage and selection for the INGARCH(p, q) model. Furthermore, the temporal aggregation and systematic sampling, which were widely studied in continuous-valued time series, have received the attention of Su and Zhu [51] in integer-valued time series.
For clarity, we have sorted out the main relationships of the models reviewed in this section in Figure 2.

Methodologies
There are many theory-oriented advances. The proofs for the existence and uniqueness of stationary distributions in the above-mentioned literature or in the earlier INGARCH literature are based on approximation techniques, the weak dependence or the theoretical framework using the Feller property, e-chain and Lyapunov's method. In contrast, to prove the existence and uniqueness of a stationary distribution and absolute regularity for nonlinear GARCH and INGARCH models under a mild assumption, Doukhan and Neumann [17] treated Z t = (X t , · · · , X t−p+1 , λ t , · · · , λ t−q+1 ) as a time-homogeneous Markov chain where {λ t } is the accompanying process of random intensities, and compensated for missing Feller properties with coupling results. Specially, besides a geometric drift condition, only a semi-contractive condition was imposed, which means a subgeometric, rather than the more usual geometric, decay rate of the mixing coefficients. This result not only enriches the theoretical proof technique, but also broadens the application area of the INGARCH model.  Further, supposed that {X t } follow from Poisson INGARCH, Neumann [52] use the contraction property twice: first, under the contraction condition of the intensity process, Neumann [52] obtained the contraction property of the Markov kernel connected with Z t in terms of a suitable Wasserstein metric, and then the existence and uniqueness of a stationary distribution follows via the Banach fixed point theorem; next, the almost effortlessly absolute regularity of the count process was established by using the contraction property once more; finally, Neumann [52] constructed a coupling of the original and the bootstrap process, and proved the existence and uniqueness of a stationary version of this joint process as well as absolute regularity of the joint count processes. Notably, the last item implies that the model-based bootstrap method proposed by Neumann [52] is more general than most of the existing papers on the consistency of bootstrap. Specifically, the bootstrap process mimics the random behavior of the original counting process, rather than being limited to the plausibility of certain specific statistics. More broadly, Doukhan et al. [53] derived the absolute regularity at a geometric rate not only for stationary Poisson GARCH processes, but also for models with an explosive trend. Recently, Aknouche and Francq [54] considered the existence of a stationary and ergodic solution of a general Markov-Switching autoregressive conditional mean model, of which the INGARCH model is one of the variants.
The contribution of Douc et al. [55] is establishing the necessary and sufficient conditions for the identifiability of observation-driven models including the pure INGARCH model and its numerous extensions, such as the pure INGARCH model with thresholds or exogenous covariates.
Recent advances in estimation methods regarding the INGARCH model are reviewed as follows. One of the main reasons for the utility of negative binomial models is that Poisson INGARCH is less flexible than models based on overdispersed conditional distributions when modeling overdispersed series. However, parameter estimation for a negative binomial INGARCH model is usually implemented based on Poisson quasi-maximum likelihood estimation (quasi-MLE or QMLE), where the pitfall is that Poisson-QMLE is likely to fail to achieve its full asymptotic efficiency with the presence of overdispersion. To clear this hurdle, Aknouche et al. [56] proposed two negative binomial QMLEs (NB-QMLEs), including the profile NB-QMLE calculated while arbitrarily fixing the dispersion parameter of the negative binomial likelihood, and the two-stage NB-QMLE consisting estimation for both conditional mean and dispersion parameters. Similarly, the two-stage weighted least square estimators (WLSEs) proposed by Aknouche and Francq [57] are general for time series data, and feasible for the INGARCH model. WLSEs can be implemented without fully specifying the conditional distribution or time series structure, and enjoy the same consistency properties as QMLEs when the conditional distribution is mis-specified, even if the conditional variance is mis-specified. Additionally, a data-driven strategy was identified to find asymptotically optimal WLSEs. This actually provides prerequisite support for further relaxation of the conditional distribution assumption for the count time series. For the model in Aknouche and Scotto [46], parameter estimation is conducted by using a two-stage WLSE, and Xu et al. [58] considered a saddlepoint MLE for a special case of this model.
The commonly used MLE method is highly influenced by outliers, so there are several works dedicated to establish robust estimation methods. For Poisson INGARCH models, Li et al. [59] proposed a robust M-estimator by using a new loss function inspired by the Tukey's biweight function. It is the construction of this loss function that contributes to this work. One of the disadvantages of Tukey's function is that it drops to zero so that the effect of very large outliers or leverage points is completely suppressed, which implies the possibility of multiple solutions to the estimated equations. Then, an intuitive idea is to construct a hybrid loss function that does not fully suppress the effects of very large outliers or leverage points. Li et al. [59] proposed a new loss function by twofold improvement. Along similar lines, Xiong and Zhu [60] introduced a robust estimation for the one-parameter exponential family INGARCH(1,1) models.
Moreover, Kim and Lee [61] used the minimum density power dispersion estimator as a robust estimator for INGARCH models whose conditional distribution belongs to the one-parameter exponential family. There are two advantages to this approach: the first is simplicity, i.e., it contains only a single tuning parameter that controls the trade-off between robustness and efficiency; the second is the ability to balance robustness and efficiency, providing considerable robustness while retaining high levels of efficiency as the tuning parameter approaches zero. Further, Xiong and Zhu [62] and Xiong and Zhu [63] used the Mallows' quasi-likelihood estimator and the minimum density power dispersion estimator as a robust estimator for negative binomial INGARCH models, respectively. For negative binomial INGARCH models, Elsaied and Fried [64] also developed several robust estimators including robustifications of method of moments and ML-estimation, one of which was an alternative to the robust estimator proposed by Xiong and Zhu [63].
Another drawback of MLE for INGARCH models is that numerical results are sensitive to the choice of initial values. Hence, Li and Zhu [65] proposed the mean targeting estimation, which is an analogue to variance targeting estimation used in the GARCH model. In addition, there have been some other advances in estimation. Jo and Lee [66] introduced the mean targeting QMLE based on INGARCH models, which provides a new perspective. When shifting to focus on more specific estimation problems, it is a common problem that the estimation performance of the intercept parameter is inferior to that of other parameters, either in the Poisson or negative binomial INGARCH model. Integrating the likelihood function by assuming a conditional distribution is one option to eliminate this obstacle. Hence, Pei and Zhu [67] adopted the marginal likelihood to estimate the intercept parameter in the negative binomial INGARCH model.
There are also some areas that are of interest to scholars. For example, to test the parameter variation of INGARCH, the cumulative sum (CUSUM) statistics has been the most popular method in recent years (Lee and Lee [68], Lee et al. [69], Lee [70], Vanli et al. [71], Weiß and Testik [72]), and Lee and Kim [73] reviewed a recent progress regarding the change point test for integer-valued time series models. Further, as with the techniques employed to present robust estimates, Kim and Lee [74] introduced a robust change point test based on density power divergence. Michel [75] considered the limiting distribution of the INGARCH(1,1) with α + β = 1.

Applications
The practical application of INGARCH models has developed considerably in recent years, especially in the period of COVID-19 when there is a strong demand for analysis of count time series data such as daily new infections in various countries or regions. This has also given rise to many valuable research topics. Agosto and Giudici [1] focused on COVID-19 contagion and digital finance, and presented Poisson INGARCH of the daily newly observed cases to understand the contagion dynamics of the COVID-19. In addition, the purpose of Agosto et al. [2] is to monitor COVID-19 contagion growth and came to the interesting conclusion that policy measures aimed at reducing infection are very useful when the it is at its peak and can reduce reproduction rates. Souza et al. [10] considered the tv-DINGARCH model with covariate/exogenous time series and applied to the daily number of deaths due to COVID-19 in Ireland. In contrast, Roy [12] focused specifically on the data for New York City because the epidemic status in this city lasted for a month and with the help of the ongoing blockade, it recovered significantly within about three months. Hence, it is this temporal variability in the data that drew Roy [12] to explore its trends by the time-varying Bayesian INGARCH model. Similarly, Giudici et al. [3] was also concerned about the time-varying features of COVID-19 and proposed Bayesian time-dependent Poisson autoregressive models. Additionally, Gning et al. [76] focused on COVID-19 in Senegal and China. Moreover, the dynamics of COVID-19 infectivity in Saudi Arabia were evaluated in Alzahrani [77] by using two statistical models, namely the log-linear Poisson autoregressive model and the ARIMA model. The results of this study showed that the log-linear Poisson autoregressive model had superior predictive performance. At the same time, many application-oriented works have actually proposed new models to meet the requirements. For example, Xu et al. [78] proposes a comprehensive adaptive log-linear zero-inflated generalized Poisson INGARCH to describe crime counts in Byron and Australia, and the features of this data set include autocorrelation, heteroscedasticity, overdispersion and excessive number of zero observations.
In addition to COVID-19, INGARCH-type models are employed in other areas such as stock trading, co-tracking of commodity marketsand so on (Chen and Khamthong [79], Agosto and Raffinetti [80], Jamaludin et al. [81], Algieri and Leccadito [82], Aknouche et al. [83], Berentsen et al. [84], Cerqueti et al. [85]). It is worth mentioning that the INGARCH-type model has been favored for human influenza research even before the outbreak of COVID-19. Specifically, Chen et al. [86] involves the INGARCH-type model when examining the causal relationship between environmental fine particulate matter and human influenza in Taiwan. The study on traffic forecasting of Kim [87] affirms the value of the INGARCH model for applications. The prediction model of Kim [87] was generated by estimating the parameters of the INGARCH process and predicting the Poisson parameters of the future step ahead process using conditional MLE methods and prediction procedures, respectively. They came to the conclusion: "INGARCH captures the characteristics of network traffic better than other statistical models, it is more tractable than neural networks (NN), overcomes the black-box nature of NN, and some statistical models perform comparable or even better than NN, especially when there is insufficient data to apply deep NN". Another application that tends to be humanistic and social Anavatan and Kayacan [88], the aim of which was to reveal the relationship between the number of femicide, female unemployment rate, male unemployment rate and the amount of information in Turkey by using INGARCH model.
In summary, the application scenarios of the INGARCH-type models can be as specific as a vehicle prediction at an intersection or as macro as a humanistic exploration of a country, and what is more valuable is that the exploration about their application still remains expected and promising.

Z-Valued Time Series
The previous section is concerned with count time series (i.e., N -valued time series). However, time series that allow both non-negative and negative integer values (i.e., Zvalued time series) are also worth investigating, whose research value is reflected in the following aspects: first, there are many typical application areas for Z-valued time series, e.g., describing score gaps in sports, trading changes in finance (Xu and Zhu [20]); secondly, the non-stationary property embodied in such time series is also of interest, such as differenced series that are initially non-stationary (Gonçalves and Mendes-Lopes [89]). Let {Z t } t∈Z be a Z-valued process, and z 1 , · · · , z T be a finite set of observations. A recent study on the mixing properties of Z-valued GARCH processes can be found in Doukhan et al. [90]. Below, we review contributions on Z-valued time series in recent years.
Firstly, the Skellam distribution is introduced into the INGARCH-type models as a conditional distribution. Alomani et al. [91] proposed the Skellam GARCH(1, 1) model defined in Definition 10, where the Skellam is the distribution of the difference of two independent Poisson variates and thus allows for both non-negative and negative integer-valued variables. The specific definition of the Skellam distribution is placed in Table 2. The benefits of the Skellam INGARCH model are that they allow time-varying variance and nonstationarity in the mean for time series, and the conditional maximum likelihood and conditional least squares methods have been developed for estimation of the parameters. Table 2. Summary of basic properties of distributions related to Z-valued INGARCH.
The drawback of the symmetric Skellam INGARCH(1,1) in Definition 10 is that only the simplest case is considered, i.e., the conditional expectation value is equal to zero. Cui et al. [92] proposes an asymmetric Skellam INGARCH to eliminate this limitation. Furthermore, the modified Skellam (MS) distribution (see Table 2) is introduced and thus the modified Skellam INGARCH in Definition 11. Specifically, Cui et al. [92] added a new parameter, γ, to the standard Skellam distribution whose mission is to control the distance between the probabilities P(Z t = 0) and min{P(Z t = 1), P(Z t = −1)} by its magnitude and sign. Hence, the modified Skellam INGARCH can flexibly compensate for the over-or under-representation of specific integers (−1, 0, 1).
Another path to modeling the Z-valued time series is to extend the N -valued IN-GARCH models by introducing a sequence of i.i.d. binary random variables independent of {Z t }, {Q t }, taking values at 1 and −1 with equal probability 0.5, such as a two-sided Poisson distribution. For example, Hu [18] and Hu and Andrews [19] proposed a Poisson Z-valued Glosten-Jagannathan-Runkle GARCH (PZG) model as follows.
Definition 12. We call {Z t } an integer-valued asymmetric GARCH process of orders p and q, if for all t ∈ Z, Hence, {X t , t ∈ Z } is a non-negative integer-valued stochastic process; conditioned on past information up to and including time t − 1, X t has Poisson distribution with mean λ t . So that {η t } and {λ t } are positive, we assume parameter α 0 > 0 and parameters α i , β j , for i = 1, · · · , p, j = 1, · · · , q, are all non-negative, with model orders p ≥ 1, q ≥ 0.
It can be seen that a two-sided Poisson distribution is employed and the structure of η t is inspired by the Glosten-Jagannathan-Runkle GARCH (GJR-GARCH, Glosten et al. [93]). The main role of GJR-GARCH here is to portray asymmetric responses in the volatility of Z-valued time series, such as the presence of leverage effects in financial time series, and is thus the highlight of the PZG model. Furthermore, Xu and Zhu [20] proposed the geometric Z-valued GJR-GARCH model based on the shifted geometric distribution that is more flexible than Poisson distribution. Additionally, Xu and Zhu [20] expanded from a two-side equality probability of 0.5 to a broader form where the ratio of positive, negative and zero values can be controlled by a parameter ρ (Definition 13). Definition 13. The geometric Z-valued GJR-GARCH model is defined as where SGe(·) denotes the shifted geometric distribution, ω > 0,|γ i | ≤ 1, α i ≥ 0, β j ≥ 0, for i = 1, · · · , p, j = 1, · · · , q, p ≥ 1, q ≥ 0. Specially, {Q * t } is a sequence of i.i.d. random variables taking values at 1, 0 and −1 with probabilities ρ, 1 − 2ρ and ρ, respectively, where 0 < ρ < 1 2 .

Bounded Count Time Series
In what follows, {B t } t∈Z consists of bounded counts with range {0, · · · , n} for given n ∈ N . This type of data also includes categorical time series, with binary time series being a special case. In terms of the INCARCH-type models, the study of {B t } differs from that of unbounded counts {X t } in two ways: one is that the candidates for the conditional distribution are required to be discrete distributions with bounded range of values; the other is that the dynamic structure of the conditional mean needs to be adjusted accordingly to the constraints of parameters of the conditional distribution, which leads to the possibility that the dynamic structure cannot be directly attached to the conditional mean, but some kind of functional transformation of the conditional mean. Accordingly, researchers have been bursting at the seams with recent innovative work in these two areas.
The definition of the bound INGARCH (BINGARCH) model is obtained by the distributional assumption that B t is generated by a bounded-count distribution and the normalized conditional mean P t = 1 n E(X t |X t−1 , · · · ) , where with the additional constraint a 0 + ∑ p i=1 a i + ∑ q j=1 b j < 1. First, for convenience, we assume that B t is generated by the conditional binomial distribution Bin(n, P t ).
Then, as with the constraints on the INGARCH-type model with linear dynamic structure mentioned previously, the BINGARCH-type model also requires constraints a 0 > 0 and a 1 , · · · , a p , b 1 , · · · , b q ≥ 0 to ensure positive P t in (15). In line with the idea of Weiß et al. [11], Weiß and Jahn [22] solved this puzzle in the BINGARCH-type model with the help of the soft-clipping functions, which enables the migration and application of this framework proposed by Weiß et al. [11]. Specially, Weiß and Jahn [22] considered the normalized conditional mean as follows: where f (·) was set to be the soft-clipping function  (16), (17) and (18). The model would be well-defined without any further restrictions, but some reasonable constraints such as |α i |, |β j | < 1 for i = 1, · · · , p and j = 1, · · · , q, and α 0 ∈ (0, 1 + p + q) are needed.
The soft-clipping binomial INGARCH was employed in Weiß and Testik [72] once more. The fascinating question considered by Weiß and Testik [72] is how the performance of the control chart is affected if the CUSUM control chart is designed based on the assumption of a completely linear data generation process, while the true one is only approximately linear. Weiß and Testik [72] aptly exploits the fact that the soft-clipping binomial INGARCH is an approximately linear counterpart of BINGARCH, and draw the conclusion that, in general, chart designs are robust to model mis-specification when parameters are specified, whereas the opposite result is obtained when the parameters are estimated.
The binomial distribution is a traditional choice for studying bounded count time series, as in (16), due to its simple form and relatively well-established properties. The motivation for the innovation against the conditional distribution is that there is a fixed relationship between the variance and the mean of the binomial distribution, denoted as the binomial index of dispersion (BID), similar to the equidispersion property of Poisson. The BID for a random variable X taking values in N is defined as BID = nVar(X) E(X)(n − E(X)) .
Additionally, the BID of the binomial random variable is calculated to be 1, which indicates that (16) is not competent for modeling data with BID > 1. Hence, Chen et al. [23] proposed a new class of INGARCH models with beta-binomial (BB) variation, which is a generalization of Chen et al. [94], in Definition 15, where the BID of BB distribution takes values in the interval (1, n). The specific form of the beta-binomial distribution is available in Table 3. To analyze the high volatility in time series counts, covariates were further introduced by Chen et al. [24], and thus a covariate-driven beta-binomial INGARCH model was proposed in Definition 16.
The beta-binomial distribution employed in Chen et al. [23] and Chen et al. [24] allow to model bounded data with under-dispersion. Then, Chen [95] turned its attention to a rare case, i.e., an under-diversified pseudo-binomial data set. It is the discrete beta (DB) distribution that competently models bounded data with under-dispersion, equivdispersion, and over-dispersion. Motivated by this and the soft-clipping function used in Weiß et al. [11], Chen [95] proposed a new soft-clipping discrete beta GARCH model as follows.
Definition 17. The soft-clipping discrete beta GARCH(1,1) model is defined by where the definition of DB is placed in Table 3, sc c (·) is defined in (18), |α 1 | + |β 1 | < 1, n bot = 0 or 1 and n top ∈ N is the predetermined upper limit of the range. Table 3. Summary of basic properties of distributions related to BINGARCH.

BB (Chen et al. [23])
The PMF of BB is The beta-binomial distribution approximately reduces to the usual bi- and The PMF of DB is n top ∈ N is the predetermined upper limit of the range and n bot = 0 or 1 is the predetermined lower limit of the range (Z taking values in {n bot , n bot + 1, n bot + 2, ..., n top }) Categorical time series are also a type of bounded count-valued time series. In the study of INGARCH-type models, categorical time series are usually presented in vector form, which can be modeled by an autoregressive multinomial logistic time series model with a latent process and is defined by a GARCH-type recursive equation. Suppose that we observe a process with state space {0, 1, · · · , n} and define a (n − 1)-dimensional vector Y t = (Y 1t , Y 2t , · · · , Y (n−1)t ) , for 1 ≤ t ≤ n, such that Y kt = 1, if the kth category is observed at time t, for all k = 1, 2, · · · , n − 1. Moreover, is defined as a vector of conditional probabilities and p t ≡ (p 1t , p 2t , · · · , p (n−1)t ) . For the last category n, set Y nt = 1 − ∑ n−1 k=1 Y kt and p nt = 1 − ∑ n−1 k=1 p kt . The dynamic structure is dependent on p t . For instance, the following linear dynamic structure is a classic: where d is a vector and A, B are matrices of appropriate dimensions. It is easy to see that the disadvantage of this modeling approach in application is the problem of using multidimensional methods to deal with univariate data resulting in increased pressure on parameter estimation and other troubles. However, its theoretical development is being gradually refined. Fokianos and Truquet [21] employed a useful coupling technique to study the ergodicity of infinite-order finite-state stochastic processes, making significant improvements to previous conditions on the stationarity and ergodicity of these models.
In addition to dealing with categorical time series in vector form, Liu et al. [96] proposed a simple and less computationally stressful way of modeling, which was essentially an innovative approach to conditional distributions from an application perspective. Liu et al. [96] first introduced a new zero-one-inflated bounded Poisson (ZOBP) distribution defined as where p 1 ≥ 0 and p 2 ≥ 0 are the inflated parameters for the states 0 and 1, respectively, with the constraint p 1 + p 2 < 1, λ > 0 is the intensity parameter and the integer M ≥ 2 is a given upper bound. This distribution is suitable for depicting data for air quality classes that are predominantly excellent, where 0 and 1 represent excellent and good air quality, respectively. Liu et al. [96] defined a new INGARCH-type model based on the ZOBP distribution. For ease of presentation, we denote the ZOBP distribution in (23) with (p 1 , p 2 ) = (0, 0) by P * (λ, M). Let {D t } be an i.i.d. sequence with the following probability distribution: where p 1 ≥ 0, p 2 ≥ 0 and p 1 + p 2 ≤ 1. Then, the ZOBP autoregressive (ZOBPAR) model is defined as follows: M) and where d > 0, a ≥ 0, b > 0, and D t satisfying (12) is independent of B * t .
Compared with the model in vector form, the ZOBPAR model is concise both in terms of estimation and its own form. The follow-up Liu et al. [97] is also based on this and is an application-oriented study completed using Bayesian estimation methods.

Multivariate Integer-Valued Time Series
It can be seen that univariate INGARCH models have been well-studied in the literature, but progress in multivariate INGARCH models has lagged somewhat in comparison. It is encouraging to note that there have been some recent developments.
We start with a review of some studies based on further improvements of the bivariate Poisson (BP) INGARCH model. The definition of the BP-INGARCH model is given here.

BP ( Cui and Zhu [27])
A bivariate Poisson distribution is defined as a product of Poisson marginals with a multiplicative factor δ, whose PMF is given by
For the BP-INGARCH model, Lee et al. [25] showed the asymptotic normality of the conditional MLE and introduced the CUSUM test for parameter change based on the estimates and residuals. Additionally, Kim et al. [26] focused on a robust estimation method for BP-INGARCH models using the minimum density power divergence estimator.
The limitation of the BP-INGARCH model is that it can only handle positive crosscorrelation between two components and is not competent for cross-correlations. A new BP-INGARCH model was proposed by Cui and Zhu [27] allowing for negative cross-correlation. Cui and Zhu [27] enabled an alternative definition of the BP distribution, just denoted as BP with definition given in Table 4, i.e., the product of Poisson marginals and a multiplicative factor δ that can promote positive, zero, or negative cross-correlation.
Further, a class of flexible BP-INGARCH(1,1) model was introduced by Cui et al. [98]. This class of models cover three distributions determined by different special multiplicative factors, making the portrayal of dependence more flexible.
For correlated bivariate count time series data, it is worth mentioning that a new flexible bivariate conditional Poisson INGARCH model was recently proposed by Piancastelli et al. [99] to capture negative and positive cross-correlations as well. Although all of these models are flexible in terms of contemporaneous correlation, the explicit form of the correlation structure of Piancastelli et al. [99] is easier to assess. Piancastelli et al. [99] also provided a detailed comparison of these methods for bivariate count time series for reference.
The next review is no longer limited to bivariate INGARCH models, but multivariate INGARCH models. Multivariate count time series remains an area of research with a vast scope, both in terms of theoretical approaches and application-oriented research. Much of the existing methodological literature does not focus only on INGARCH models, but is more broadly applicable to various types of count time series models; see Fokianos [100] for some recent methodological developments including multivariate INGARCH models. Moreover, Fokianos et al. [101] introduced an overview of statistical analysis for some models for multivariate discrete-valued time series based on higher-order Markov chains, where several extensions are highlighted including non-stationarity, network autoregressions, conditional non-linear autoregressive models, robust estimation, random fields and spatiotemporal models. We only show the contribution made by the latest literature Lee et al. [102] because it was not included in Fokianos [100]. Definition 21. Let Y t = (Y t1 , · · · , Y tm ) , t ≥ 1, be the time series of counts taking values in N m , and p i (y|η) = exp{ηy − A i (η)}h i (y), y ∈ N , which stands for the probability mass function of one-parameter exponential family, wherein η is the natural parameter, A i (η) and h i (y) are known functions, and both A i and B i = A i stands for the derivative of A i , are strictly increasing. We then consider the following model.
where B i (η ti ) = M ti , f θ is a non-negative function defined on [0, ∞) m × N m depending on the parameter θ ∈ Θ ⊂ R d for some d = 1, 2, · · · , and η t = (η t1 , · · · , η tm ) : Unlike other authors who have devoted more effort to specifying the joint distribution of multivariate time series and the marginal distributions of their components, Lee et al. [102] argued that the conditional mean equation forms the bulk of the modeling and that the specification of the underlying joint distribution is not a major concern. We think this premise is a reasonable presupposition. It is well-known that the INGARCHtype model is built based on two types of assumptions, i.e., assumptions of conditional distribution and dynamic structure. Then it is intuitively reasonable that the contemporaneous dependence of the multivariate count time series can be reflected in the INGARCH model in two ways: the joint distribution and the multivariate dynamic structure. This presupposition reduces the complexity of modeling. Specially, although each component of Y t is modeled using a univariate INGARCH model in Definition 21, the dependence structure is imposed by the conditional mean process M t .
Another emerging area is the development of INGARCH-type models applicable to time-varying network data, which can be considered as a special class of multivariate count time series. Therefore, a part of recent research has been devoted to establish INGARCH models and their methodologies applicable to time-varying discrete network data. Let {Y t = (Y 1t , · · · , Y Nt )} be a network with N nodes. The structure of the network is completely described by the adjacency matrix A = (a ij ) ∈ R N×N , where a ii = 1 for any i = 1, · · · , N and a ij = 1 that means the presence of a directed edge from i to j, a ij = 0 otherwise. In general, any time-varying discrete network data whose relationship can be modeled by an adjacency matrix can be considered as a multivariate time series.
In a similar vein to the development of univariate integer-valued time series, the Poisson distribution and linearity assumptions continue to be used instinctively in order to establish a universal framework. Armillotta and Fokianos [28] considered the Poisson network autoregressive (PNAR) models for count data with a non-random adjacency matrix. We show here only the simplest linear form of order 1: where n i = ∑ i =j a ij is the out-degree, i.e., the total number of nodes, which i has an edge with.
The PNAR(1) model reduces the inference complexity by incorporating network information into the dependence structure, where the response of each individual can be explained by its lagged values and the average effect of its neighbors in (27). Note that Equation (27) does not include information about the joint dependence structure of the PNAR(1) model. It is then convenient to rewrite (27) in vector form, where β 0 = β 0 1 N ∈ R N , with 1 = (1, 1, · · · , 1) ∈ R N and the matrix G = β 1 W + β 2 I N , where W = diag{n −1 1 , · · · , n −1 N }A is the row-normalized adjacency matrix, with A = (a ij ), so w i = (a ij /n i , j = 1, · · · , N) ∈ R N is the i-th row vector of the matrix W, satisfying W ∞ = 1, and I N is the N × N identity matrix. {N t } is a sequence of independent N-variate copula-Poisson processes. The main methodological contribution of Armillotta and Fokianos [28] was the study of the asymptotic properties of such models by employing L p -near epoch dependence and α-mixing, rather than based on the assumption that i.i.d. on which the development of all network time series models discussed so far has strongly relied. Further, Armillotta et al. [29] reviewed some of the work by Armillotta and Fokianos [28] and provided a unified framework for the statistical analysis of both continuous and integer-valued data with a known adjacency matrix. Armillotta and Fokianos [28] also specified a log-linear PNAR model for the count processes, and another recent work Armillotta and Fokianos [30] was closely related to this, where a quasi-score linearity test for continuous and count network autoregressive models was developed.
It can be seen in (27) and (28) that the PNAR model assumes that all individuals are homogeneous and they share a common autoregressive coefficient. This is a somewhat detached assumption from reality. Therefore, Tao et al. [31] proposed the grouped PNAR, which divides individuals into different groups and describes heterogeneous node behavior with group-specific parameters. Compared to the original PNAR model, the constraints are relaxed while competently portraying the heterogeneity. Specially, all individuals could be classified into K groups in the setting of Tao et al. [31], and each group was characterized by a specific set of positive parameters θ k = (ω k , α k , ρ k , β k ) ∈ R 4 , for 1 ≤ k ≤ K.

Definition 23.
A grouped PNAR model can be constructed as for each i = 1, · · · , N, and t ≥ 1. Following the PNAR model, the parameters ω k , α k , ρ k , β k represent the group-specific baseline effect, regression coefficient on past observations, network effect, and regression coefficient on past intensity processes, respectively. Note that we assume the adjacency matrix A is asymmetric, which covers the special case of symmetric networks. To distinguish between groups, latent variable z ik ∈ {0, 1} was defined for each object i, where z ik = 1 if object i is from the k-th group, and z ik = 0 otherwise. Assume {(z i1 , · · · , z iK ) , 1 ≤ i ≤ N} is a sequence of i.i.d. multinomial random vectors with number of events n = 1 and probability γ = (γ 1 , · · · , γ K ) . Here, γ k represents the group proportion satisfying γ k ≥ 0 and ∑ K k=1 γ k = 1.
Tao et al. [31] explored the accuracy of model estimation and prediction when the group labels were unknown and the number of group K is mis-specified, respectively. There is already a body of mature research on network data, but it is still a relatively emerging topic in the INGARCH field. Therefore, many further attempts to consider time-varying networks from the perspective of count time series are worthwhile, such as optimization of conditional distributions, nonlinear dynamic structures, and related hypothesis testing, or time-varying network data with upper bounds on the number of edges.
In addition to modeling and methodology, INGARCH is popular for applicationoriented analysis of time-varying networks. For example, Agosto and Ahelegbey [103] used a financial network model to study the contagion effects between business sectors based on discrete data, and tested the conditional means (and volatilities) of default counts across economic sectors estimated by Poisson INGARCH and their dependence in shocks. Through an empirical analysis of corporate defaults in Italy over the period 1996-2018, a high degree of intersectoral vulnerability was concluded by Agosto and Ahelegbey [103], in particular at the onset of the global financial crisis in 2008 and in subsequent years. Such a wide range of application prospects is accompanied by a desire for theoretical development.

Discussion and Conclusions
The purpose of this section is to present some potentially useful research topics based on the methodology and applications reviewed in the previous sections.
(1). First, we focus on the softplus function sp c (·) and soft-clipping function sc c n (·), which contribute to the modeling of unbounded counts and bounded counts allowing for negative auto-correlation, respectively. For the sake of clarity, sc c n (·) is used next as an example. As already mentioned, the advantage of sc c n (·) is that the support set is R and is nearly linear on (0, n], which allows the parameters in the dynamic structure of the BINGARCH model not to be restricted to positive numbers. This is certainly an excellent innovation, and one that seems worth exploring further. The images of sc c n (x) corresponding to different parameters c or n are reported in Figures 3 and 4. It is obvious that the slope of sc c n (x) is small when x < 0 and tends to zero as x decreases. Moreover, the parameter c has a small moderating effect on this tendency, while n even has almost no effect. This insensitivity to negative values may lead to concerns that the corresponding INGARCH models do not fairly model positively and negatively correlated data.  The reasons for this confusion can be summarized as follows. Although negative parameters are allowed to appear in the dynamic structure of the conditional expectation, such leniency seems to be somewhat offset by the fact that P t (we use Bin(n, p t ) as an example here) that were obtained from sc c n (x) are concentrated around 0 when x < 0, while sc c n (x) is approximately linear when x > 0. The sc c n (x) does allow for the existence of negative correlation, but the extent of the negative correlation is to be explored further. However, to put it another way, it seems to us that sc c n (x) provides a new perspective of modeling zero-inflated counts. In contrast to the previous idea of embodying the zeroinflated feature in a conditional distribution, it seems possible to assign this task to a dynamic structure through sc c n (x). (2). There is still a large demand for innovations dedicated to bounded count time series. We personally think that there are two feasible directions for exploration: one is the innovation of conditional distribution, and the other is to continue to deepen the research based on vector forms. From an application point of view, it is worth exploring to truncate existing distributions in a reasonable way, or to find distributions that themselves take values in bounded sets of integers. For the study of vector forms, where the theory is relatively well-developed, it is imperative to ease the estimation pressure.
(3). The emergence of new methodologies or research topics in other fields or in the broader field can stimulate the development of the INGARCH-type models. For example, Pedersen and Rahbek [104] presented the theory for testing for reduction of GARCHX type models with an exogenous covariate to standard GARCH type models. Many INGARCHtype models, including some of the recent literature mentioned earlier, are also focusing on covariates. Thus, the tests and methodologies proposed by Pedersen and Rahbek [104] can actually inspire existing INGARCH models to achieve tests of the reasonableness of introducing covariates. Similarly, refer also to Aknouche and Francq [105] and Debaly and Truquet [106].